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P^PL^SMA  THEORY  AND  SIMULATION 
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A.  Lower  Hybrid  Drift  Instability  2d  Simulations  \  \  ^ 

■  \  Yu-Jiuan  Chen  (Prof.  C.K.  Birdsall,  Dr.  W.  Nevins,  LLL) 

We  have  improved  our  method  of  loading  particles  in  x,y  so  that  the  initial  density  n(x,y) 
Is  as  near  our  theoretical  equilibrium  as  possible. /  As  a  result,  only  a  small  amount  of  electron 
heating  was  observed  in  the  first  few  cycles  of  u^per  hybrid  oscillations. 

We  have  used  small  drift  velocities  (  <  v„  )  in  our  simulations  to  date,  allowing  elec¬ 

trostatic  fields,  zero  beta.  Unfortunately,  this  choice  makes  both  the  growth  rates  and  satura¬ 
tion  levels  of  the  lower  hybrid  drift  instability  small,  hence  very  difficult  to  simulate.  We  are 
using  the  multibeaming  quiet  start  loader  for  ion  particles.  Both  multibeam  instabilities  and  the 
lower  hybrid  drift  instability  have  been  observed. 

A  new  diagnostic,  gyrokinetic  phase  scattering  plots'  have  been  added  to  study  the  elec¬ 
tron  heating.  We  have  been  successfully  subtracting  the  ExB  drift  corresponding  to  all  the  low 
frequency  waves  from  the  electron  velocities.  In  our  simulations,  a  small  amplitude,  slowly 
growing  upper  hybrid  wave  was  observed.  We  found  the  electron  gyro-kinetic  energy,  including 
electron  slashing  energy  in  the  upper  hybrid  wave,  increasing  in  time.  However,  comparing 
with  ion  temperature,  electrons  are  still  very  coldUv 

More  study  of  nonlinear  mechanisms  will  oe  continued.  Also  simulations  with  larger  drift 
velocities  will  be  tried,  so  as  to  emphasize  the  effects  sought. 


Magnetized  Multi-Ring  Instabilities 


Niels  Otani,  Dr.  B.I.  Cohen  (LLL),  and 
Dr.  M.J.  Gerver  (MIT)  (Prof.  C.K.  Birdsall) 


I.  Continuous  rings  (M.J.  Gerver) 

We  wish  to  find  an  approximate  simple  expression  for  the  stability  threshold  of  magnet¬ 
ized  plasma  whose  distribution  function  consists  of  N  monoenergetic  components,  arranged  to 
approximate  a  Maxwellian.  The  dispersion  relation  is  given  by  Equation  (5)  on  p.  11  of  QPR 

II,  1980,  but  with  a  factor  of  l/N  in  front  of  the  ion  term,  which  was  inadvertently  left  out  of 
Equation  (5).  (The  solution  of  Equation  (5)  is  correct  in  Figure  5,  p.  13  of  QPR  II,  1980; 
checked  by  N.  Otani.]  Near  w  -  / ft,  (and  this  is  qualitatively  valid  as  long  as  |w  —  /ft,|  <  ft,) 
we  can  neglect  all  but  one  term  in  the  sum  over  /,  and  the  dispersion  relation  becomes, 

to  -  /ft,  ,  _  Upti  ft,2  _  2  <£  «*<*  ,  ,  kLv, 

/ft,  1  +  ~  k "  n- 

'  « 

2  4  ft,2  .  2*jr.| 

~  (1) 


where  we  have  used  the  large  argument  forms  of  J\  and  Jt  to  obtain  the  last  expression.  The 
arguments  of  the  sine  function  for  two  consecutive  values  of  s  differ  by 


Empirically,  the  instability  first  appears  at  fcjvyjVfl,  =  2,  which  means  that  the  different  terms 
in  the  sum  over  s  are  out  of  phase  by  ~  2  radians,  and  the  sum  over  all  s  is  typically  about 
equal  to  one  of  the  terms  in  the  sum. 

£  n,2  2klVs  _  n }  __  i 

n,  J~  klv,2  "  4N2  (2) 


(Note  that  this  result  depends  on  regular  spacing  of  the  rings.  If  the  values  of  v,  were  chosen 
randomly,  with  the  correct  distribution  function,  then  the  sum  would  be  greater  by  V77). 


Instability  occurs,  roughly,  when 


a >  —  Id  j  ^  Cl  j 

~  / Q,  J_ 
/ft,  Z  l 


v-ja^  j_ 

in,  *  2  n 

This  last  inequality  is  based  on  the  empirical  observation  that  instability  first  appears  at 
/  ~  IN  12. 

Combining  (1)  and  (2)  yields 

<■>  ~  'to/  f.  |  A,2  __  1 


<■>  -  /n, 

/n,  i  +  ne2  a>2 

O-lCli  __  U>1  [ 

“*■  ~ ; - n?|1  + 

Combining  (3)  and  (4)  gives  the  condition  for  stability: 


This  result  differs  by  N  from  the  rule  given  in  QPR  III,  p.  2,  1980;  some  reconciliation  is 
needed.  In  practice,  we  must  use  an  "effective  N’,  which  is  the  number  of  terms  in  the  sum 
over  s  which  contribute  significantly  to  the  sum.  This  number  will  be  quite  a  bit  less  (maybe 
by  a  factor  of  6)  than  the  total  number  of  terms  N.  In  principle,  however,  otjJCl}  should  go  as 
N1  for  sufficiently  large  N.  For  N  sufficiently  large  [probably  something  like  3 (m,/m,)w  in 
practice],  there  will  be  no  instability. 


2.  An  Electrostatic  Dispersion  Relation  for  a  Uniform  Magnetized  Plasma  Having  a  Discrete- 
Particle  or  Rings-and-Spokes  Velocity  Distribution  (Niels  F.  Otani  and  Bruce  I.  Cohen) 

We  continue  to  examine  the  electrostatic  instability  associated  with  particle  velocity  distri¬ 
butions  found  in  plasma  simulation  computer  codes.  In  the  previous  QPR  the  analysis  of  this 
instability  was  based  on  the  behavior  of  the  equilibrium  distribution  /0  -  Xan8(vj.-vin).  We 

n 

now  extend  the  analysis  to  include  the  discreteness  of  particles  in  gyrophase  space  ("rings  and 
spokes").  As  before,  the  equilibrium  is  assumed  to  be  spatially  uniform  and  uniformly  magnet- 
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The  derivation  in  (a)  is  due  to  Otani.  Cohen  offers  an  alternative  derivation  in  (b)  based 
on  action-angle  variables  that  confirms  the  analysis  in  (a)  but  is  simpler  and  more  compact. 
Conclusions  are  presented  in  (c). 

(a)  Derivation  of  the  Dispersion  Relation 

We  choose  as  our  phase-space  coordinates  p,  9',  X,  Y,  z,  and  v.  where,  in  terms  of  the 
usual  cartesian  phase-space  coordinates, 


Vx  +  Vj? 

2n 

(1) 

y 

9'  —  arctan  —  —  Sit 

(2) 

vy 

X~x+n 

(3) 

V  X 

Y~y~H 


eBn 


Here  ft  — - .  Note  the  explicit  time  dependence  of  9'  on  t. 

me 

We  wish  to  consider  equilibria  of  the  form 

/o  ■  /oO*,vr,0') 


(4) 


(5) 


These  are  valid  equilibria  since  p,  v.,  and  9'  are  all  constants  of  the  motion.  In  particular,  9'  is 
the  gyrophase  measured  relative  to  a  corotating  reference  angle  and  hence  is  constant  along  an 
equilibrium  trajectory. 

The  linear  dispersion  relation  is  obtained  by  considering  the  perturbed  distribution 

/  -  /oO*,vz,0')  +  8f(p,vz,9\X,  Y,t )  (6) 


Notice  in  this  derivation  we  are  only  considering  modes  with  spatial  dependence  in  the  plane 
perpendicular  to  the  magnetic  field  B0 z.  Also  it  is  assumed  these  modes  are  collisionless: 


4L 

dt 


-0 


(7) 


This  equation  combined  with  the  equations  of  motion 


r  —  -£-8E  +  — vxBo 
m  me 


yields  the  linearized  equation 


+  —rvi*8Ei  +  — tt-SExy^z-1  —■  -  0 


(8) 


(9) 


dt  mil  L  L  mil  1  2fi  d9' 

Here  the  assumptions  8E  -  -V<£  (electrostatic  assumption)  and  B/Bz  -  0  have  been  used. 
Fourier  transforming  the  time  dependence  we  have: 


/cnS/(z,<u)  — 


mil 


%tS  *  + -k 


2 p  B9' 


(10) 


in  which,  for  brevity,  z  -  (p,vz,9',X,  Y).  It  is  interesting  to  note  that  the  explicit  time  depen¬ 
dence  of  Vj_(z,/)  leads  to  the  coupling  of  different  frequency  components: 


8/(z,»)  - 


2mu 

-  e-'f  S£-(z ,w-n) 


,tu+fl) 

a/o 

.  _L 

2m 

a/o  _ 

/ 

a/o 

Bp 

2m 

B9‘ 

a/o 


(id 

(12) 


a 


Here  rt  =  vJSl  and  &£*  —  5£x  ±  i8Ey.  Another  equation  for  8E  is  obtained  as  follows 


8E(z,«>)  -  -  fdt  e‘“'(V<*)(z ,/) 


where 


(Vd»)(z ,/)  -  Vd>(x,y,t)  x  -  x  +  ,tc<*(0+n,) 

y  -  Y  -  yin  W'+fti) 

Defining  xL  -  (x,y),  k^  -  (kx',ky'),  and  X  -  (X,Y),  we  have 


n  C  d*ki  pda)1  /  <k .  '** .  -•/) ,  i\ 

V«(Xi,r)  -J  -^p-J  —  e  --  /Mfki.oO 


Now  substituting  back  into  Eq.  (14): 

*(i,r)  -  J*-^  e'ki‘xe'kirt (16) 

where  fcx'  -  — lk/1  sin  a'  and  /c/  -  — |kj.'|  cos  a'.  With  the  help  of  the  usual  Bessel  identity, 
we  obtain 

6E(z,a)  -  -  J* e'kl  Xy£tJr(k/ rg) e,ri*'~a) ikL<j>(ki ,u>+l' Cl)  (17) 

A  third  equation  may  be  found  from 

p(x,t)  -  fdz  en<fi3(x  -  r(z,r))8/(z,r)  (18) 

where  n0  is  the  unperturbed  number  density,  p  is  the  perturbed  charge  density,  and  r(z,r)  are 
the  spatial  cartesian  coordinates  as  functions  of  z  and  t.  Fourier  transforming, 

p(k,«)  -  en0J dt  dz  (19) 

Using  Poisson’s  equation  and  the  Bessel  identity  again, 

*(k,a»)  -  -^^2. / dz  e~*>2 e''ki'XXy, (kLrt)  (Trf,#'-“)8/ (z,a>-/ ft )  (20) 


where  A:x  —  — |kj.|  sin  a  and  -  — |kj  cos  a.  Substituting  (12)  and  (17)  into  (20)  we  obtain 

m, 

_<(/'-/-  DW-n)  Q  if-  ;  3*- 

+ 2 —  v  -  2i-^rk^+(r-,-i)ft)|  (23) 

where  mp  =  (4tt n0e2/ m)'h  and  gQ  =  J* dv.  f0. 

The  determinant  of  the  implied  (infinite)  matrix  serves  as  the  dispersion  function  for  the 
general  case  /o  -  fo(jJ-,vz,9').  Equation  (23)  assumes  only  one  contributing  species  of  particles; 
generalization  to  more  species  is  straightforward.  Notice  that  the  eigenmodes  for  this  system 
do  not  in  general  have  the  simple  time  dependence  e~,a“,  but  instead  are  composed  of  a  mix¬ 
ture  of  frequencies  ft  apart.  This  is  a  consequence  of  the  time-dependent  equilibrium  /0.  As  a 
check  we  find  that  the  well-known  gyrophase-independent  dispersion  relation 

1  +  £f  2*  -  0  (24) 

klJ  dp  i  w  -  /ft 

is  recovered  when  dg^dS'  -  0.  Also  Eq.  (23)  is  found  to  be  consistent  with  the  reality 


condition  <t>  (kifa»)  —  <£(— kj,,-<i»). 

Now  choose  g0(fi,9')  to  be  the  discrete-particle  velocity  distribution 

*o(M')  -  IMO*  -  M.)  1 8(«*  -  (rjp-  +  9  „)) 


m—  l 


where  a„,9„,  fi and  M„  are  constants.  In  this  case  Eq.  (23)  reduces  to 


Kl  I.Hj  * 


JjMA9~a) 


kj„ 


M/. 


( l+jM„)Ji'Ji+jM  +  UiJ'i+jM 
x - i_ - — (kj.,<u— fl ) 

<i>  —  in 


(25) 


(26) 

(27) 


The  argument  of  all  Bessel  functions  is  £]/„;  r„  is  the  gyroradius  corresponding  to  fi„. 

For  the  rings-and-spokes  distribution,  M„  =  M  and  9n  =  0  for  all  n.  Equation  (27) 
may  now  be  written  as 


where 


/ 


AJJ+*J 


T  D(l+iMAJ) 

ft2  /  ^L_y 

n 


(28) 


(29) 


and 


2>(/,Ay)  -  +  (l+MAj)Jl+MljJ,']  (30) 

We  find  that  A  is  hermitian  when  <n  is  real  and  see  that  only  terms  such  that 
M|8y|  <  2A:i(rB)m„  contribute  since  Jt(kLrn)  goes  to  zero  very  rapidly  once  |/|  is  appreciably 
larger  than  kjH.  In  particular,  if  M  >  2ki(r„)mtx  only  diagonal  terms  contribute.  In  other 
words,  the  rings-and-spokes  distribution  behaves  essentially  like  the  corresponding  gyrophase- 
independent  distribution  when  the  wavelength  is  longer  than  the  spacing  of  M  particles 
arranged  on  a  circle  of  radius  (/■*)„«• 

(b)  Alternative  Derivation 

An  alternative  and  somewhat  simpler  derivation  is  based  on  the  use  of  action-angle  vari¬ 
ables  as  described  in  Aamodt  et  al.  (1981).  The  Vlasov  equation  can  be  written  as 

fit)  -  /(0)  -  j*  dt'  \HJ),  (31) 


where  H  is  the  Hamiltonian  and  ( }  is  the  Poisson  bracket.  The  single  particle  Hamiltonian  for 
a  flute  mode  (k’B  -  0)  is 


H  -  Hq  +  H\  —  nCl  +  qi  exp  (He c  -  ioit)  +  c.c.  - 

fin  +  qdt'£J„(kvJd)  exp  (in 9  -  io>t)  +  c.c., 

n 


(32) 

(33) 


where  X  and  Y  are  the  x  and  y  guiding  center  positions,  fi  =  mvf/2fl,  9  -  arctan  vx/vy, 
O  =  qBo/mc,  B0  is  a  uniform  external  magnetic  field  in  the  z  direction,  and  k  has  been  taken 
parallel  to  x  without  loss  of  generality.  The  Vlasov  equation  can  be  linearized  with  respect  to 
Hx/Hq: 


•  fa//, 


e/o 


a//i  e/o  1 


/  *5  A  \ 


t 
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In  Eq.  (34),  the  equilibrium  distribution  function  f0  has  been  taken  to  be  independent  of  X 
and  Y ,  as  we  assume  the  plasma  to  be  uniform. 

Canonical  transformation  to  the  new  variables  (0,/Z)  using  the  generating  function 
Fi  *»  (0  —  gives  0  —  bF^biJi  —  9  —  {It,  —  /I  —  bF^b9,  and  H  —  H  +  bF^Jbt  —  e<t> , 
where 

e<f>  -  «<££./,,  exp  [//r(0+n/)  —  not]  +  c.c.  (35) 

n 

To  zero  order  in  the  perturbation  e4>>  the  equations  of  motion  are  trivial, 

-£(9,]I-,X,Y)(0)  -  0;  (36) 

all  of  the  independent  variables  are  now  constants  of  the  zero  order  motion.  A  suitable  equili¬ 
brium  distribution  function  composed  of  Nt  angular  spokes  is  given  by 

fo  -  noNf1  £8(0  -  9j)f0(bar).  (37) 

/- 1 

Equation  (34)  becomes 

,  a/o  .  .9J„  a/o 

nJ„—zr  +  i—z — =" 

fx  exp  ( ikx  —  /to/)  —  —  q4>  £  - — - ^ exp  (//i(0+n/)  —  /to/].  (38) 

,  to  —  nil 

Use  of  Poisson’s  equation  — V20  —  4 ir^qsfd3v  fls  enables  us  to  obtain  a  dispersion 

1 

relation 

k2d>  -  4*r£  q,  J dJLd9  fXs  - 


4w£/»o02/  d?f  o(m>Z~ ^ 

J  t,n 


"TrW  -  (n-/)/,-^- 

_ OU 


(O  —  flfi 


*  #«■'£  exp  [/(rt-/)(0;+n/)]^ 

j 

If  the  angles  9j  are  evenly  spaced,  0y  -  (J-l)2n/Nf,  then 

(Vr'X  exP  I/(«-/)Wy+ftrtl  -  exp [/(«—/) fl/]8„t/+/vtfll, 

j 

where  8,y  is  the  Kronecker  delta  and  m  —  0,  ±  1,  •  •  • .  Hence, 

k2d>  -  k2£0,  exp  Or  fl/)  -  £  47r«0q2J’  /o(m>£  4>p 

r  s  p 

n^pU,J„)  ~  ( n~l)J ,^r 
x  I~  ^  to-pn-/,n  ^  »-./+■*,  e»p 


(39) 

(40) 

(41) 

(42) 

(43) 

(44) 


from  which  follows 

4Trn0?2 


*,-£ 


k2 


J  dfifoG) 


i  y  _ 

'  ”  (<o—rCl)—nfl 


n~~T  ~  ™NiJ„-my-—r 

+  I  ^ 


«.« 

m#0 


to  —  r(l  —  n  Cl  +  mNpCl 


(45) 


(46) 
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Equation  (28)  in  (a)  is  easily  recovered  from  this  expression. 

The  discrete  evenly  spaced  angles  have  produced  an  alias  coupling  of  temporal  modes 
separated  in  frequency  by  mN9Sl,  where  m  -  0,±1,±2,  •  •  •  As  argued  in  (a),  we  can  arrange 
for  the  coupling  to  be  weak  and  thus  recover  the  desired  dispersion  relation  for  an  equilibrium 
with  no  gyrophase  structure, 

4  4  "IF'® 

'  -  -p- 2 ^ /•© I  ■  W> 

if  the  coefficients  multiplying  &r-mNt  for  m  0  are  negligible,  so  that  only  the  coupling  of  4>r 
to  itself  survives  in  Eq.  (46).  For  small  argument,  the  Bessel  function  has  the  limiting  form 
Jv(z)  —  (z/2)p/i>\,  and  \Jv(z) I  «  1  for  jz|  <  v/2.  In  a  typical  simulation  with  a  thermal 
velocity  distribution,  max  (vj  <  4v,A.  Hence  the  argument  of  the  Bessel  functions  in  Eq. 
(46)  satifies  the  inequality  z  —  kvJCl  <  (4A:maxvrt/fl).  The  coupling  of  <f>r  to  will  be 

negligible  if 


N»  > 


(48) 


The  analysis  of  Kim  and  Otani  in  QPR  II,  1980  describing  the  multi-beaming  instability  of  a 
warm  magnetized  plasma  made  up  of  many  rings  forming  a  Maxwellian  velocity  distribution  can 
then  be  applied  to  Eq.  (47).  Generally  speaking,  for  NVj_  equally  weighted  rings,  stability  can 
be  achieved  for 


*i>i 


OJr 


(49) 


Thus,  a  quiet-start  simulation  ought  to  be  stable  to  magnetized  multi-ring  instabilities  if  the 
inequalities  in  Eqs.  (48)  and  (49)  are  both  satisfied. 

(c)  Conclusions 

In  conclusion,  dispersion  functions  for  the  general  distribution  fo(n,v2,0'),  the  discrete- 
particle  distribution,  and  the  rings-and-spokes  distribution  have  been  derived  and  found  to  be 
consistent  with  well-known  gyrophase-independent  dispersion  relations.  We  find  that  the  eigen- 
modes  are  in  general  not  simple  sinusoidal  functions  of  time  but  are  instead  a  mixture  of  fre¬ 
quencies.  However,  in  the  case  of  the  rings-and-spokes  distribution  it  was  shown  that  the 
different  frequencies  decouple  when  the  number  of  particles  per  ring  (i.e.,  number  of  spokes) 
becomes  somewhat  larger  than  Mr,)™*  in  which  case  the  dispersion  relation  reduces  to  its 
gyrophase-independent  counterpart.  With  a  sufficient  number  of  discrete  9'  and  ju.  classes 
(spokes  and  rings)  it  should  then  be  possible  to  perform  quiet-start  simulations  that  are  stable 
to  magnetized  multi-ring  instability. 

(d)  Acknowledgments 

We  are  grateful  to  W.  M.  Nevins  for  a  number  of  illuminating  discussions  and  to  C.  K. 
Birdsall  for  his  interest  and  encouragement.  The  research  of  B.  I.  Cohen  was  performed  under 
the  auspices  of  the  U.  S.  Department  of  Energy  at  Lawrence  Livermore  National  Laboratory 
under  Contract  W-7405-ENG-48.  Reference 

R.  E.  Aamodt,  B.  I.  Cohen,  Y.  C.  Lee,  C.  S.  Liu,  D.  R.  Nicholson,  and  M.  N.  Rosenbluth, 
Phys.  Fluids,  24,  55,  (1981). 


C.  Long  Field- Reversing  Ion  Layers 


Douglas  S.  Hamed  (Prof.  C.K.  Birdsall) 
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(1)  Paper 

The  paper  following  (3)  was  presented  at  the  Third  Symposium  on  the  Physics  and  Tech¬ 
nology  of  Compact  Toroids  in  the  Magnetic  Fusion  Energy  Program  held  at  Los  Alamos, 
New  Mexico,  Dec.  2-4,  1980.  It  summarizes  work  on  ion  layer  kink  instabilities  done 
through  the  end  of  October  1980. 

(2)  Code  Development 

The  following  modifications  have  been  made  during  November  and  December  to  our 
two-dimensional  hybrid  code,  AQUARIUS:  (a)  The  mover  has  been  vectorized,  produc¬ 
ing  a  factor  of  2.7  increase  in  efficiency,  (b)  Finite  electron  pressure  has  been  included, 
(c)  A  new  version  has  been  developed  which  employs  conducting  wall  boundaries  and  has 
the  capability  to  handle  vacuum  regions. 

(3)  Results 

Two  new  results  have  been  obtained  from  AQUARIUS:  (a)  New  diagnostics  have  enabled 
measurement  of  the  real  part  of  the  frequencies  of  ion  layer  kink  instabilities.  These  fre¬ 
quencies  have  been  found  to  agree  with  theoretical  predictions,  (b)  Saturation  of  ion 
layer  kink  instabilities  has  been  found  to  be  due  to  heating  of  the  layer.  The  heating, 
which  occurs  during  instability  growth,  causes  the  layer  to  thicken.  As  the  layer  becomes 
thicker  the  betatron  frequency  of  the  layer  particles  is  reduced.  Once  the  betatron  fre¬ 
quency  has  become  sufficiently  low,  then  the  self-magnetic  field  index  drops  below  the 
instability  threshold  value  and  the  instability  ceases  to  grow. 


Kink  Motion  of  Long  Field-Reversing  Ion  Layers 

Douglas  S.  Haraed 

Electronics  Research  Laboratory 
University  of  California 
Berkeley,  California  94720 

Kink  instabilities  have  been  studied  in  long  field-reversing  ion  layers.  The  configuration 
is  shown  in  Figure.  1.  An  ion  beam,  at  radius  R  and  of  thickness  a ,  crosses  an  external  mag¬ 
netic  field  The  beam  current,  Jsb9,  produces  a  self-magnetic  field,  5/,  which  tends  to 

reverse  the  total  magnetic  field  on  axis.  The  layer  is  immersed  in  a  uniform  background  plasma 
such  that  nb«np,  where  nb  and  np  are  the  densities  of  the  beam  and  plasma,  respectively. 
There  is  assumed  to  be  no  variation  in  the  axial  direction  (i.e.,  9/9z— 0).  Kink  modes  in  these 
layers  correspond  to  perturbations  of  the  form  e*-«,exp  ( in0—ia>t)r ,  with  n^2. 

Long  layer  kink  modes  have  been  studied  by  Lovelace1.  Analytic  results  were  obtained 
for  uniform  layers  with  sharp  boundaries  within  the  approximations  (a/R):<<l. 
W/ift)2«l,  and  vi=constant,  where  Cl  is  the  layer  rotation  frequency  and  vA  the  plasma 
Alfven  velocity.  Stability  thresholds  were  obtained  for  two  cases:  (1)  co2«(vA/a)2  and  (2) 
w2<(vA/a)2  with  w2«a>2,  where  u,  and  «,  are  the  real  and  imaginary  parts  of  the  frequency. 
In  both  cases  the  condition  r)s<n2—\  was  found  to  be  sufficient  for  stability,  where  rjj  is  the 
self-magnetic  field  index.  If  {  is  defined  to  be  the  loading  factor  (f=|2?/(r—0)/2?/|),  then  -qs  is 
of  order  (2 R/a)C 

In  order  to  have  a  useful  comparison  for  our  simulation  results  it  was  necessary  to  gen¬ 
eralize  the  preceding  theoretical  analysis  by  eliminating  the  assumptions  a f<(vA/a)2,  w;<o>,2, 
and  vj=constant.  In  addition  we  chose  to  study  the  stability  of  a  rigid  rotor  distribution.  The 
approximations  were  eliminated  by  performing  a  numerical  solution  of  Lovelace's  equations. 
With  the  less  severe  constraints  on  frequencies  our  numerical  results  produced  more  restrictive 
requirements  on  the  self-magnetic  field  index  for  stability.  Minimum  values  of  r)s  for 
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instability  are  shown  in  Table  1.  It  can  be  seen  that  the  thresholds  increase  with  increasing 
background  density. 

A  two-dimensional  simulation  code,  AQUARIUS,  has  been  developed  and  applied  to  ion 
layer  kink  motion.  The  code  treats  ions  as  particles  and  electrons  as  a  massless  fluid.  The 
Darwin  version  of  Maxwell’s  equations  are  used  (i.e.  the  transverse  displacement  current  is 
neglected).  Maxwell’s  equations  are  coupled  with  the  electron  momentum  equation  (with  iner¬ 
tial  terms  neglected)  by  the  assumption  of  quasineutrality  to  determine  the  time  evolution  of 
the  electric  and  magnetic  fields.  The  ion  particles  obey  the  equations  of  motion.  Standard 
particle- in-cell  techniques  are  employed  to  determine  forces  and  accumulate  currents. 

The  code  is  initialized  with  an  equilibrium  exponential  rigid  rotor  ion  layer  distribution. 
Instabilities  have  been  observed  corresponding  to  n -2,3,4,  and  5  modes.  Figure  2  shows  the 
particle  positions  in  their  initial  state  for  a  layer  with  5=1.4.  In  Fig.  3  the  particle  positions  of 
this  same  layer  have  become  significantly  perturbed  after  an  n— 3  mode  has  grown  to  large 
amplitude.  When  linear  growth  rates  from  simulations  are  compared  with  the  predictions  of 
our  numerical  extension  of  the  theory  in  Ref.  1,  excellent  agreement  is  obtained.  The  simula¬ 
tion  results  and  the  corresponding  theoretical  growth  rates  for  an  ion  layer  with  R/a= 5  are 
shown  in  Figs.  4  and  5.  The  maximum  growth  rate  for  a  given  mode  occurs  at  n=a)p/(l, 
where  <op  is  the  layer  betatron  frequency.  Figure  5  shows  how  the  growth  decreases  with 
increasing  background  plasma  density.  Nonlinear  effects  have  been  found  to  halt  the  linear 
growth  of  kink  instabilities.  The  saturation  amplitudes  decrease  with  increasing  azimutha'  mode 
number.  The  instabilities  have  been  observed  to  grow  and  saturate  without  destroying  the 
field-reversed  layer.  Kink  instabilities  instead  result  in  a  thicker  field-reversed  state  with  many 
non-axis-encircling  particles.  The  final  state  appears  to  be  stable  but  with  higher  electric  fields 
and  substantial  electron  currents. 


1.  R.V.  Lovelace,  Phys.  Fluids  22,  708  (1979). 
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TBL.  1.  Instability  threshold  values  of  the 
self-magnetic  field  index. 
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FIG.  2.  Initial  particle  positions  in  the  /— 0 
plane  for  an  ion  layer  with  {—1.32. 


FIG.  3.  Particle  positions  after  ten  ion- 
cyclotron  periods.  An  it—  3  instability  has 
grown  to  a  large  amplitude  and  nonlinear 
effects  have  become  important. 
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FIG.  4.  Growth  rates,  w„  for  ion  layers 
with  (rt/a)= 5  and  (v^/cJ-.OOHS  as  a 
function  of  the  loading  factor. 
(Q  n— 2.A  »-3,n«-4,0"-5). 


FIG.  5.  Growth  rales  of  «- 2  modes  for  ion 
layers  with  (R/a)~ 5  and  {—.65  as  a  func¬ 
tion  of  background  plasma  Alfven  speed. 


Section  II 

CODE  DEVELOPMENT  AND  MAINTENANCE 


A.  Implicit  Particle  Simulation  using  Moments  with  Orbit  Averaging 


V.  Thomas  (Dr.  B.  I.  Cohen  LLL  and  Prof.  C.  K.  Birdsail) 

A  first  attempt  at  including  orbit  averaging  in  electrostatic  implicit  moment  particle  simu¬ 
lations  has  failed.  The  motivation  for  attempting  to  include  orbit  averaging  in  a  moment  impli¬ 
cit  code  is  for  the  improvement  of  particle  statistics.  (  See  the  previous  QPRs  for  implicit 
moment  simulation  references  and  orbit  averaging  references.  )  The  algorithm  is  shown  in  Fig¬ 
ure  1.  Starting  at  time  NkT  with  £'v+l*  we  advance  the  partictes  from  N\T  to  (jV+1)A7\ 
Then  we  use  the  current  density  averaged  over  the  micro  steps  from  (N— 1/2) AT  to 
(N+1/2)&T  and  the  charge  density  averaged  over  the  microsteps  from  NIT  to  {N+D&T  to 
predict  the  charge  density  at  time  (JV+3/2)AT  via 

ps+m  -  (p}*+1'2  -  ;-v+1A  T  (50) 


and 


r 


,V+1 


l/J 


s  _  Ha 
ma 


bx 


{pj*+I/J-qJ/.J"+,/2£>r 


(51) 


the  brackets  indicate  orbit  averaging  and  the  tilde  indicates  a  predicted  quantity.  P  is  the  pres¬ 
sure,  alpha  is  the  species  label,  £  is  a  linear  combination  of  EN  ,  £*+l,  and  £'v+2  ,  and  n 
represents  the  particle  density.  The  field  at  time  (Af+2)AT  is  solved  for  implicitly  by 


—(!—«)  V2<**+2  -  eV2<£*+I  -  pN+3/2  (52) 


where  p'v+3/2  is  predicted  from  the  first  two  equations,  e  is  a  biasing  parameter  intended  to 
control  dissipation. 

This  algorithm  is  unstable  ,  even  for  small  time  steps.  We  are  working  to  understand  the 
results  and  we  are  persuing  alternatives. 


/^Vion 


NAT 


(N  +  DAT 


'V, 


electron 


l 

(N+2)AT 


Figure  1.  Using  Es*[  the  ions  and  the  electrons  are  advanced  on  possibly  different  time  steps  up  to  the  time 
(Af+DAT.  Then  £V+J  is  solved  for  implicitly  and  the  cycle  repeats. 
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B.  A  NEW  METHOD  FOR  IMPLICIT  PARTICLE  SIMULATION 

Alex  Friedman*.  A.  Bruce  Langdon**,  and  Bruce  I.  Cohen** 

The  following  is  a  brief  account  of  our  recent  work  on  implicit 
methods  for  particle  simulation.  The  algorithm  we  describe,  named  BAAL 
(for  Bay  Area  ALgorithm),  was  developed  during  the  final  quarter  of 
I960,  while  one  of  us  (A.F.)  was  a  postdoctoral  associate  with  the 
U.  C.  Berkeley  Plasma  Theory  and  Simulation  Group.  This  work,  which  has 
application  to  bath  inertial  and  magnetic  confinement  fusion  plasma 
simulation,  is  continuing  at  the  Lawrence  Livermore  National  Laboratory, 
with  continuing  close  ties  to  the  U.  C.  Berkeley  group. 

The  first  section  below  (written  primarily  by  A.B.L. )  has  been 
excerpted  from  our  contribution  to  the  1980  LLNL  Laser  Fusion  Annual 
Report.  It  provides  an  overview  of  our  work,  and  describes  an  ideal 
(gridless)  form  of  the  algorithm.  The  second  section  below  briefly 
describes  the  cloud-in-cell  algor ithra  used  in  the  testbed  code  developed 
at  U.  C.  Berkeley.  The  final  section  presents  some  illustrative 
results.  Details  of  this  work,  other  forms  of  the  BAAL  algorithm,  and 
our  analyses  of  the  time-differencing  algorithms  will  be  presented  in 
future  publications. 

*Present  affiliation:  Lawrence  Livermore  National  Laboratory 

**Af f i 1 iation:  Lawrence  Livermore  National  Laboratory 


1.  Ideal  Algorithm 


We  have  made  several  advances  in  the  quest  for  more  efficient 
simulation  of  low-frequency  plasma  phenomena.  Our  most  adaptable  and 
reliable  tools  for  study  of  kinetic  plasma  behavior  are  the  "particle" 
codes,  but  the  stability  of  these  codes  has  previously  required 
resolution  of  the  electron  plasma  period  in  the  time  integration,  even 
when  the  phenomenon  under  study  took  place  on  the  much-longer  ion  time 
scale . 

Analogous  limitations  on  time  step  in  other  problems,  such  as  heat 
flow  and  chemical  rate  equations,  are  overcome  through  the  use  of 
implicit  time  integration  schemes.  In  particle  codes,  although  implicit 
methods  have  been  analyzed  theoretically,1  their  application  has  been 
inhibited  by  the  very  large  number  of  nonlinear  equations  to  be  solved 
simultaneously,  about  equal  to  the  number  of  particle  coordinates  plus 
the  number  of  zone  quantities  (electric  and  magnetic  fields). 

We  have  begun  to  experiment  with  a  new  method  for  solution  of  these 
equations  in  two  steps.  In  this  method,  equations  are  first  set  up  for 
the  fields.  Although  these  equations  involve  information  from  the 
particles,  the  number  of  equations  to  solve  simultaneously  is  only  the 
number  of  field  quantities,  defined  on  the  zones.  Since  the  number  of 
zones  is  much  less  than  the  number  of  particles,  and  the  resulting 
matrix  equations  are  sparse  and  well  conditioned,  their  solution  is 
convenient  with  known  methods.  Once  the  fields  are  known,  the  particle 
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coordinates  can  be  readily  solved  for  serially,  one  particle  at  a  time. 
Here,  we  will  discuss  simulations  having  only  the  electrostatic  field. 

First  we  consider  f inite-dif ferenced  equations  of  motion  for  the 
particles  which  have  the  necessary  stability  at  large  time-step  and  are 
accurate  for  the  low  frequency  phenomena  to  be  studied.  In  Ref.  3,  a 
centered  expression  for  the  time  derivative  is  employed.  Use  of  this 
expression  leads  to  a  numerically  stable  algorithm  at  large  wpA t,  but  by 
itself  does  not  damp  high  frequency  components  of  the  motion,  which  are 
aliased  into  the  Nyquist  frequency  -  stable  odd-even  oscillations  of 
large  amplitude  can  persist.  In  Ref.  2,  the  derivatives  are  biased 
using  a  scheme  such  as 

v»+rvn  *  <  f*n+l*4an-l  >At 

*n+l~*u  *  <  K-l+K-l  >At 

This  scheme  damps  unwanted  high-frequency  oscillations,  while  low- 
frequencies  are  very  weakly  damped,  as  desired:  ( Imw)/«=©(wA t )3 .  Such 
schemes  are  members  of  a  class  whose  application  to  plasma  simulation 
has  been  analyzed  in  detail  in  Ref.  1  and  in  our  present  work.  We  have 
also  devised  a  different  class  of  schemes,  whose  simplest  member  is 

2(  va+l/2~ yn-l/2  )  "  (  va-l/2~va-3/2  )  “  *n+lAt 
xb-m'-xii  *  ya+l/2^ 


(la) 

(lb) 


(2a) 

(2b) 
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This  scheme  has  the  same  order  of  accuracy  as  Eq.  (l)>  while  requiring 
less  storage  of  time  levels.  The  presence  of  the  acceleration  at  only 
the  n+1  time  level  increases  the  damping  of  high-frequency  oscillations. 
The  optimum  design  of  these  difference  equations  is  the  first  issue  in 
practical  implementation  of  large  time-step  methods. 

In  all  implicit  schemes  the  new  positions  xa+1  depend  on  the 
accelerations  aa+1  due  to  the  electric  field  Ea+1 .  But  this  field  is 
not  yet  known,  as  it  depends  on  the  density  pn+1  of  particle  positions 
-{xn+1}’.  The  solution  of  these  coupled  particle  and  field  equations  is 
the  other  major  implementation  issue. 

In  the  method  developed  by  Denavit  and  Mason  for  this 
solution,  the  fields  at  the  new  time  level  are  predicted  by  solving 
coupled  field  and  fluid  equations,  in  which  the  pressure  tensor  is 
approximately  evaluated  from  particle  velocities  known  at  the  earlier 
time.  After  the  fields  are  known,  the  particles  are  advanced  to  the  new 
time-level,  and,  if  desired,  an  improved  pressure  tensor  is  calculated 
and  the  process  iterated. 

It  is  also  practicable  to  predict  the  future  electric  field  En  +  1 
quite  directly  by  means  of  a  linearization  of  the  particle-field 
equations.  Since  this  method  and  its  implementation  in  the  experimental 
code  BAAL  have  not  been  described  previously,  we  outline  the  concept 


here . 
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,  as  given  by  an 


The  position  xn+1  of  a  particle  at  time  level  tn+1 
implicit  time  integration  scheme,  can  be  written  as 

*n+l  38  0At2an  +  l  +  xn+l* 


where  and  is  the  position  obtained  from  the  equation  of 

motion  with  the  acceleration  aQ+1  neglected.  Since  x£®|  depends  only  on 
positions  and  accelerations  at  times  tn  and  earlier,  it  is  known.  In 
its  simplest  form,  the  BAAL  algorithm  is  derived  by  linearization  of  the 
particle  positions  relative  to  x£®] . 

One  can  regard  the  actual  position,  xa+1,  as  xa+!  P^us  a 
displacement  dx  =  /3Atzan+1.  We  form  a  charge  density  from 

th®  actual  charge  distribution  is  then  p£°|,  plus  the  change  dp 
brought  about  by  displacing  particles  by  the  amount  dx  =  xn+1-x£°]. 
Linearized,  this  increment  to  p  is1 


dp  =  -V-[p<°,)(x)dx(x)]  (4) 

To  the  same  order  of  approximation,  the  displacement  dx(x)  of  all 
particles  with  x£°]  a  x  is  obtained  with  aa+1  evaluated  at  x,  i.e. 

dx(x)  a  /SAt2(q/m)En+,(x)  (5) 


We  then  have 


dp(x)  =  -V'£x(x)EJJ^.1(x)] 


(6) 


19 

where  the  effective  susceptibility  is 

X(x)  =  /SCp^+!(x)q/m]At2  =  0<y2(x)At2  (7) 

Note  that  x  depends  only  on  the  particle  positions  ■C*n2}}‘.  and  not  at 
all  on  velocity  information  as  required  in  the  moment  equation  methods. 

With  these  two  source  contributions,  the  Poisson  equation  becomes, 
in  rationalized  c.g.s.  units, 

*’Ea+i  »  pJSl  -  *-<xEa+1)  (8) 

or 

-V[l+x]Vpn+1  =  pi®!  (9) 

This  elliptic  equation  is  solved  by  standard  methods.  The  field  -V<Pn+\ 
and  (5)  are  then  used  to  calculate  the  positions  {xn+1}.  This  algorithm 
is  reminiscent  of  the  method  for  solution  for  the  vector  potential  A  in 
some  magneto— induct ive  plasma  simulation  codes.4,5 

Should  it  be  necessary,  we  have  shown  how  to  refine  the 
approximations  used  above  by:  linearization  about  a  more  accurate 
prediction  of  *n+l  than  x£®|t  iteration,  and  a  more  accurate  evaluation 


In  our  spatial-difference  representation  of  these  equations,  Eq. 

(4)  becomes  the  gradient  of  a  zonal  p  with  respect  to  particle  position, 
and  6x  in  Eq.  (5)  depends  on  the  electric  field  in  two  zones,  using  the 
usual  interpolation.  In  this  way  we  are  assured  that  the  density  pn+1 
of  final  particle  positions  -(xn+1}  satisfies  the  code's  representation 
of  the  Poisson  equation,  ~^2Pn+1  -  pa+1»  s0  that  desirable  features 
built  into  the  time-differencing  scheme  will  be  realized  in  practice. 

We  are  delineating  the  limitations  of  implicit  methods  and  learning 
how  to  surmount  them.  For  example,  it  is  observed  that  an  unphysical 
cooling  of  the  electrons  occurs  when  uncentered  equations  of  motion  are 
used.  Our  analysis  provides  guidance  to  design  algorithms  which  retain 
desirable  dissipation  of  high-frequency  oscillations  while  minimizing 
this  unwanted  cooling.  This  side  effect  is  related  to  the  damping  of 
law  frequency  oscillations;  both  are  smaller  in  our  schemes  above  with 
third-order  damping  than  in  lower-order  algorithms. 

To  date,  the  BAAL  approach  has  been  verified  in  application  to  ion 
acoustic  oscillation  and  two-stream  instability,  with  wpAt  in  the  range 
2.5  —  4.0.  Numerical  stability  and  correct  dispersion  of  the  testbed 
code  have  been  verified  for  a  cold  simulation  plasma  up  to  cj  At  =  25. 
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II.  Cloud-in-Ce 1 1  Algorithm 

The  testbed  computer  simulation  program  we  have  developed,  QAAL, 
employs  well-known  particle-in-cell  techniques  to  assure  a  smooth 
interpolation  of  the  zonal  electric  field  onto  the  particle  locations 
and  of  the  particle  charge  onto  the  charge-density  grid.  We  illustrate 
the  scheme  with  a  ful ly-impl ic i t  (backward  differenced)  model  algorithm; 
the  actual  code  is  far  more  general. 

In  the  model  algorithm,  each  particle  k  is  advanced  via: 

v*+1  =  v“  -+  (q/rn)AtEm+1(x“+1) 
xj+l  -  x“  +  AtvB+1. 

For  each  cell  j,  we  write  the  charge  density  as: 

=  (q/Ax)IS(x®+1— Xj ) ,  (11) 

where  the  "shape  function"  S  is  the  finite-size  particle  analogue  of  the 
Dirac  delta  function,  and  the  sum  is  taken  over  all  particles  k  v.'hich 
overlap  cell  j  to  any  degree.  The  calculation  of  Em+1  requires 
knowledge  of  pB+1  before  we  have  explicit  knowledge  of  the  x“+1.  In 
essence,  the  BAAL  algorithm  calculates  this  charge  density  implicitly  by 
approximating  i t  as  a  linear  function  of  Em+1.  We  compute  a  simple 
approximation  xk  (this  is  x£°]  of  the  previous  section)  to  x“+1  via: 


(10a) 

(10b) 
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Xk  =  x®  +  AtvJ?,  (12) 


which  for  this  simple  model  is  just  a  "free-streaming"  location.  We 
expand  the  shape  function  in  a  Taylor  series: 

S(x£+1-Xj)  *  S^-Xj)  +  (xg+1-xk)as(xlc-xj  )/axk .  (13) 

and  note  that  the  expansion  is  exact  for  linear-spline  particles  •which 
do  not  cross  cell  boundaries  as  a  result  of  the  new  acceleration, 
i.e.  which  intersect  the  same  cells  at  x®+1  and  xk-  We  approximate  the 
coefficient  of  the  first-derivative  term  via 


x£+1-ik  =  (q/nOA^E""-1^1) 

*  (q/m)At2IS(Jk-xi)E“+l,  (14) 

with  a  relative  error  of  order  wTrM>AtS  *  <iEAt2/mL3cale  **  l’  where 
Lacale  is  a  typical  scale  length  over  which  E  varies.  This 
approximation  imposes  a  limitation  (shared  by  the  moment-equation 
method)  on  Ljemle  and  |E|  for  the  simulation  to  be  valid.  We  are 
considering  modifications  to  the  algorithm  which  may  eventually  allow  us 
to  relax  this  restriction. 

Using  Eqns.  (11).  (13),  and  (14),  the  field  equation  solved  is: 


(V‘E)j  ~  ^j(conv)  + 


(15) 
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where  the  "conventional"  charge  density  is: 

^J(coxit)  =  (16) 

and  the  "weights"  Wt  j  couple  each  cell  to  its  neighbors  when  particles 
are  present : 

Wjj  *  (q2At3/m)p(xk~xi)aS(xk-xj  )/3xk.  (17) 

For  linear  splines  Wj  ^=0  whenever  1 i — j | > 1 .  Since  for  linear  splines 
dS/3xk  =  +1/Ax  (or  zero),  there  is  considerable  redundancy  in  the 
calculation  of  the  Pj(conT)  and  the  Wj ^ ,  and  it  turns  out  that  (at  least 
for  the  unmagnetized  Id  electrostatic  code)  the  number  of  operations 
needed  to  weight  all  necessary  quantities  to  the  grid  is  the  same  as  the 
number  needed  in  the  conventional  ESI  code.  Since  the  field 
interpolation  and  charge  deposition  steps  tend  to  dominate  the  particle 
processing  time  on  a  vector  computer  like  the  Cray-1,  the  algorithm  is 
potentially  very  efficient. 

At  present,  the  field  equation  is  written  in  terms  of  the 
electrostatic  potential  <p ,  and  is  solved  by  successive  overrelaxation. 
Direct  solution  is  more  difficult  to  generalize  to  higher  dimensional  problems. 
Spatial  filtering  is  provided  for  through  use  of  a  fast  Fourier  transform  on  the 
potential  obtained  after  the  field  equation  has  been  solved. 


2.4 


III.  Code  Results 

As  a  first  check  of  code  performance,  we  initialized  a  cold 
electron  plasma  (with  immobile  ions)  in  a  run  with  parameters  cjp=1.0, 
At=10.0  (run  "CLD2").  For  this  run,  the  BAAL  testbed  used  a 
time-centered  version  of  the  model  difference  scheme  described  in  this 
section.  Thirty-two  cells,  and  128  particles,  were  employed,  and  there 
was  an  initial  excitation  of  mode  1  (the  longest  wavelength  mode). 

Figure  (l)  shows  a  "history"  of  the  energy  in  this  mode  as  a  function  of 
time.  An  oscillation  with  a  period  of  160  is  clearly  evident  in  the 
figure.  That  is  not  evident  is  that  the  oscillation  is  in  fact  a 
modulated  odd-even  oscillation  in  time,  i.e.  with  a  frequency  differing 
from  the  Nyquist  frequency  ff/At  by  only  the  small  shift  <5&;=27t/160  .  To 
understand  this,  we  observe  that  the  dispersion  relation  for  this 
"trapezoidal  rule"  scheme  is: 

2/WpAt  =  cot  wAt/2  (18) 

—  tan  (n/2  -  wAt/2), 


Thus,  the  mode  frequency 
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u  %  ff/At  -  4/wpAt2,  (19) 

so  that  the  frequency  shift  is  =  4/o>pAt2  =  0.04  leading  to  a  period 
of  157.  Similar  behavior  was  observed  for  At=25.0. 

When  the  difference  scheme  of  Eqn.  (1)  is  employed,  the  odd-even 
oscillations  (aliased  plasma  wave)  are  observed  to  damp  out  rapidly, 
with  70b8^*~9 -O'7  •  This  compares  well  with  the  predicted  damping  factor 
^pr«<Ts“^ •  Similar  damping  is  obtained  when  "stiffly  stable"  schemes 
such  as  that  of  Eqn.  (2)  are  employed. 

As  a  first  model  problem,  we  consider  the  electron— electron 
two-stream  instability  arising  from  the  relative  motion  of  two  cold 
beams  of  infinite  spatial  extent.  The  first  run  made,  "EE00",  was  a 
benchmark,  using  the  conventional  ESI  leapfrog  simulation  algorithm  and 
a  small  timestep,  At=0.25  with  a>p(  total  )=1 .0  .  For  this  run,  a  grid  of 
32  cells  was  employed,  and  128  particles  were  used  to  represent  each 
beam.  No  spatial  filtering  was  employed  (the  capability  had  not  yet 
been  implemented),  and  thus  the  fastest  growing  mode  was  that  with  the 
shortest  wavelength,  mode  number  16.  Figure  (2)  is  a  plot  of  the  field 
energy  as  a  function  of  time  for  this  run;  the  analytically  predicted 
growth  rate  for  the  parameters  chosen,  normalized  to  the  plasma 
frequency,  is  7pre<j=®•'^,  wkile  the  observed  growth  rate  is  7ot,s=0-29. 
Thus,  despite  the  fact  that  the  mode  in  question  has  a  wavelength  of 
only  two  cells,  there  is  good  agreement  with  simple  theory. 


For  comparison  with  the  above  run,  another  ("EE03”)  was  made  using 
the  BAAL  testbed  code.  The  timestep  was  ten  times  as  large,  At=2.5,  and 
a  time-centered  version  of  the  model  algorithm  of  this  section  was 
employed.  Thus,  there  was  no  damping  of  the  aliased  plasma  oscillations 
at  the  Nyquist  frequency  (these  would  not  be  visible  in  a  field  energy 
plot  anyhow,  since  only  the  square  of  the  electric  field  enters  into 
that  diagnostic).  The  results  of  this  run  are  shown  in  Fig.  (3).  The 
instability  is  not  evident  until  a  later  time  in  the  BAAL  simulation; 
presumably  this  is  due  to  a  smaller  initial  level  of  excitation  of  the 
unstable  mode.  Once  visible  growth  begins,  it  is  clear  that  the  growth 
rates  in  this  run  and  the  previous  run  are  very  nearly  equal,  while  the 
saturation  level  in  the  implicit  code  run  is  only  slightly  lower  than 
that  of  the  ESI  run.  There  is  thus  good  agreement  with  both  simple 
theory  and  the  benchmark  run,  despite  the  poor  resolution  of  the  mode  in 
question  on  the  grid. 

As  a  second  model  problem8,  we  consider  an  ion-acoustic  traveling 
wave,  with  wpe=1.0,  At=3.0  (run  ”1A04'').  For  this  run,  the  time-level 
biasing  of  Eqn .  (l)  was  employed,  as  was  a  spatial  filtering  which 
discarded  all  but  modes  1-4,  and  which  filtered  modes  3  and  4  rather 
heavily.  Mode  2  was  excited,  there  were  128  grid  cells,  the  mass  ratio 
mj/me  was  16,  and  there  were  2048  electrons  and  512  ions.  The  system 
length  L  was  16n,  the  electron  thermal  velocity  vthe  was  0.25  (the  ions 
were  cold),  the  Debye  length  was  0.25,  and  the  wave  vector  k  was 
0.25.  The  sound  speed  c#  was  0,0845,  and  the  wave  frequency  w  was  0.0161 
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leading  to  a  period  v  of  390.  Figure  (4)  is  a  "snapshot"  of  the 
smoothed  ion  density  as  a  function  of  x  at  the  beginning  of  the  run. 

The  corresponding  quantity  at  time  t=210  (slightly  more  than  half  a  wave 
period  later)  is  shown  in  Fig.  (5).  One  predicts  that  the  lefthnnd  peak 
should  have  moved  210cs=13.59  units  during  this  period  of  70  timesteps, 
and  measures  a  traveled  distance  of  19.3-6.1=13.2.  Thus,  the  phase 
velocity  of  the  wave  is  essentially  correct.  At  later  times  modr  4 
comes  up  out  of  the  noise  and  the  clean  structure  visible  in  these 
snapshots  disappears. 

We  have  successfully  modeled  ion  acoustic  standing  waves  with  mass 
ratios  up  to  160C  and  timesteps  At  as  large  as  4.0. 7  Up  to  6000 
timesteps  have  been  employed  in  these  runs,  which  are  "noisier”  than  the 
one  described  above  despite  the  fact  that  only  a  single  mode  was 
al lowed. 

Our  future  plans  include  the  modeling  of  warm  two-stream 
instabilities  wherein  the  physically  fastest-growing  mode  is  well 
resolved  in  the  simulation,  and  the  verification  of  the  correct  behavior 
of  a  plasma  expanding  into  vacuum.  Practical  application  of  these 
methods  to  problems  in  inertial  and  magnetic  confinement  fusion  will  be 


forthcoming . 
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C.  POLARES 


Niels  F.  Otani  (C.  K.  Birdsall) 

The  master-slave  structure  of  the  POLARES  code  described  in  the  last  QPR  is  now  in 
place,  although  rough  spots  still  remain.  We  find  the  new  code  reproduces  results  obtained 
from  the  old  code.  The  predicted  graphics  flexibility  appears  to  have  been  realized,  though 
some  bugs  remain.  Also  difficulty  has  been  encountered  in  the  organization  of  the  history  file, 
and  corresponding  portions  of  the  code  will  probably  have  to  be  rewritten.  Restart  and  negative 
timestep  capabilities  (see  previous  QPR)  have  not  been  installed  but  the  prospect  of  their  suc¬ 
cessful  implementation  remains  promising. 

We  have  learned  from  A.  B.  Langdon  (private  communication)  that  the  master-slave 
organization  used  in  this  code  is  not  the  most  efficient  since  the  code  must  generally  wait  to 
regain  core  each  time  a  switch  between  master  and  slave  is  made.  In  this  regard  and  also  from 
the  point  of  view  of  ease  in  coding,  the  use  of  overlays  probably  would  have  been  preferable. 
However,  the  present  format,  having  already  been  written,  will  be  maintained. 


D.  Solver  Updates 

H.  Stephen  Au- Yeung 

(1)  New  Parameters 

Several  new  parameters  have  been  added  to  SOLVER  to  allow  the  user  to  override  the 
labels  and  to  override  the  lower  and  upper  bound  of  the  graph;  they  are: 

Hori  Horizontal  label. 

Vertre  Vertical  label  for  the  real  part. 

Vertai  Vertical  label  for  the  imaginary  part. 

Headre  Heading  label  for  the  real  part. 

Headai  Heading  label  for  the  imaginary  part. 

Xmin  Lower  bound  of  x. 

Xmax  Upper  bound  of  x. 

Remin  Lower  bound  of  the  real  part  of  the  root. 

Remax  Upper  bound  of  the  real  part  of  the  root. 

Aimin  Lower  bound  of  the  imaginary  part  of  the  root. 

Aimax  Upper  bound  of  the  imaginary  part  of  the  root. 

Wmin  Lower  bound  of  the  second  argument  w. 

Wmax  Upper  bound  of  the  second  argument  w. 

(2)  Get  around  the  limitation  of  1024  points  for  function  plotting 

When  plotting  a  function  vs.  one  argument,  one  could  plot  as  many  as  2048  points  by 
specifing  the  following  in  the  input  file  of  tty: 

nx  —  2048  nw—  1 

E.  Quiet  Start  Method  Comparisons; 

Random,  Bit  Reversed,  and  Fibonacci  Numbers 


H.  Stephen  Au-Yeung,  and  Yu-Jiuan  Chen  (Prof.  C.K.  Birdsall;  Dr. 
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A.B.  Langdon,  LLL) 

Quiet  starts  for  Maxwellian  distribution  lor  other  /(v)  or  n0(x)]  use  some  set  of  uniform 
numbers  for  inverting  the  cumulative  distribution  function.  The  set  may  be  regular  or  other. 
We  have  chosen  to  make  some  comparisons  among  sets  stimulated  in  part  by  the  work  of 
Denavit  and  Walsh  (1980)  using  Fibonacci  numbers  based  on  that  of  Wick  (1979,  1980). 

An  earlier  method  used  by  A.B.  Langdon  in  ESI,  subroutine  INIT  (since  June  1978)  is  to 
invert  f0(v )  using  uniform  numbers  from  slow  to  fast,  hence  highly  correlated,  and  then  to 
scramble  positions  [ na  (x)  —  constant]  in  order  to  reduce  correlations  using  a  bit  reversed 
method. 

The  main  reason  for  using  the  quiet  start  loader  is  to  reduce  the  noise  level.  The  noise 
will  be  increased  by  particle  bunching.  Therefore,  we  decided  to  compare  the  uniformity  of  the 
numbers  generated  by  the  Fibonacci  number  method,  the  bit  reversed  method,  and  computer 
(CRAY-1)  generated  random  numbers.  Keep  in  mind  that  uniformity  (lack  of  bunching)  is  not 
the  only  test,  as  correlations  that  excite  various  (<o,k)  preferentially  are  also  unwanted;  tests  of 
the  latter  kind  have  not  been  made  as  yet. 

Fibonacci  numbers  are  defined  as 


<*«+!  ~  aH  +  «»-i  .  n  *■  •  •  • 


(1) 


with  ao  —  <*1 
cles  with  N  ■ 


■  1.  The  sequence  is  1,  1,  2,  3,  5,  8,  13,  etc.  Following  Wick,  we  use  N  parti- 
a„  to  attempt  to  fill  a  unit  square  uniformly  by  using 


x. 


O'  ~  j 


1,  •  •  •  ,N-an 


mod  1 


(2) 


where  (x„y,)  is  the  location  of  the  i'h  particle.  Figure  1  has  (x,y)  scatter  plots  for  N  -  233, 
377,  610,  987,  1597,  2584,  4181,  and  6765.  There  is  no  mystery  here;  the  x,'s  step  across 
slowly  and  the  y' s  step  up  rapidly,  both  uniformly,  creating  the  2d  space  lattices  shown.  The 
resolution  of  these  patterns  required  photo  printing  (fiche);  if  done  on  the  Versatec,  with  less 
resolution,  spurious  bunching  resulted,  sort  of  Moire  patterns.  The  Fibonacci  method  clearly 
produces  highly  correlated  (xjO.  For  example,  suppose  the  particles  were  allowed  to  free 
stream  as 


x'  -*  x  —  yt  modulo  1 

y'-y 

treating  y  like  velocity.  Using  S  -  987,  note  that  at  t  =  1,  there  will  be  roughly  21  nearly  per¬ 
fect  vertical  bars,  such  that  if  (x,y)  were  (x,vx)  there  would  be  a  large  density  produced  in 
mode  21,  defeating  the  pu^ose  of  the  quiet  start.  At  other  values  of  t,  other  Fibonacci  "cry¬ 
stals"  will  "resonate"  similarly. 

For  the  bit-reversed  method,  we  set  N  «■  2"  where  n  is  an  integer.  The  x,  are  chosen  as 
before,  with  x,  -  (/  4-  Vi)  /  N,  with  i  running  0  to  N- 1.  The  y,  are  chosen  using  the  bit- 
reversed  algorithm  in  ESI  plus  a  1/2 N.  The  scatter  plots  are  shown  in  Figure  2,  for  N  -  64, 
128,  256,  512,  1024,  2048,  4096,  and  8192.  This  method  appears  to  our  eyes  to  avoid  correla¬ 
tions  in  small  (x,>)  regions  (no  obvious  lattice)  although  there  are  a  few  closely  spaced  (x- 
y's).  (There  appears  to  be  symmetry  about  ±45®.)  Hence,  the  resonance  possibility  or  magni¬ 
tude  appears  less  here. 

Bit  reversed  fractions  were  not  familiar  to  us,  so  a  short  explanation  and  algorithm  is 
given  here.  The  table  below  makes  bit  reversed  fractions  from  the  binary  equivalent  of  the 
x,’s. 


k 


Figure  2.  xny,  space  with  the  y,  chosen  from  a  bil  reversed  sequence. 
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n 

N 

Decimal 

Binary 

Binary  reversed 

Fraction 

0 

1 

xo-0 

0. 

.0 

NS 

o 

1 

O 

1 

2 

X(  -  1 

1. 

.1 

>-,  -  1/2 

x2-2 

10. 

.01 

-  1/4 

2 

4 

x3  -  3 

11. 

.11 

y3  -  3/4 

*4  —  4 

100. 

.001 

y4  -  1/8 

xs-5 

101. 

.101 

y5  -  5/8 

x4-6 

110. 

Oil 

>6  “  3/8 

3 

8 

x?  —  7 

111. 

.111 

OO 

1 

e- 

x, -  8  1000.  .0001  >-8  —  1/16 

x,-9  1001.  .1001  y9  —  9/16 

xlo-10  1010.  .0101  >-(0-5/16 

x„-ll  1011.  .1101  >-,,-13/16 

X[j  —  12  1100.  .0011  >-12-3/16 

X(3  —  13  1101.  .1011  >-,3-11/16 

xw  —  14  1110.  .0111  y,*  —  7/16 


4  16  xls  -  IS _ till. _ .1111  y,5  -  15/16 


5 

32 

16-31 

loooo.  -  inn. 

.00001  - 

.11111 

1/32  -  31/32 

6 

64 

32-63 

100000.  -  111111. 

.000001  - 

.111111 

1/64  -  63/64 

The  values  of  y  used  in  ESI  are  those  in  the  table  plus  1/2 N.  A  readable  code  for  producing 
the  bit  reversed  (yt  +  1/2 N)  follows. 

integer  Num,  n,  K,  KthBit,  TwoToK 
real  B 
namelist  /in/  n 
c 

c  Read  in  n 
c 

read(59,in)  n 
c 

c  Bit  reverse  algorithm  for  n  bits 
c 

Num  -  2**(n-l)  -  1 
do  10  I  -  1,  2**<n-l) 

Num  —  Num  +  1 
c 

c  Reverse  Num  into  B. 

c  The  "and"  function  performs  a  bit-by-bit  logical  and 
c  to  its  arguments  (Num  and  TwoToK  in  this  case), 
c 

B  -  0. 

do  20  K  -  0,  N-l 
TwoToK  -  2"*K 
KthBit  —  and(Num,TwoToK) 

B  -  2.*B  +  KthBit/TwoToK 
20  continue 

c 

c  B/(2**n)  here  is  the  I’th  number  in  the  bit  reverse  sequence 
c 

write(59,9010)  B,  B/(2**n) 
c 

10  continue 

call  exit 
c 

9010  format  ( fB.0,  4x,  fl5.8  ) 

end 
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end 


For  n  —  1,  this  produces  1  and  bit  reversed  fraction  plus  Ml S  -  1/2 

2  1,  3  and  1/4,  3/4 

3  1,  5,  3,  7  and  1/8,  5/8,  3/8,  7/8 

4  1,  9,  5,  13,  3,  11,  7,  15  and  1/16  etc. 

5  1,  17,  9,  25,  ,  31  and  1/32  etc. 

Both  the  Fibonacci  and  bit  reversed  methods  tend  to  fill  y ,  space  uniformly  with  only  the 
first  few  x,’s  (or  with  any  partial  set  of  the  x,’s).  This  is  the  point  of  the  uniformity  seen  in  the 
figures.  Hence,  if  v,’s  are  selected  sequently  (say,  from  inverting  /(v)  exp-v2/2v,2], 

meaning  highly  ordered,  then  using  these  methods  to  scramble  the  x’s  will  tend  to  produce  the 
full  /(v)  over  a  small  region  is  x.  This  filling  obviates  the  need  for  dividing  x-space  into  sub 
regions  each  filled  with  /(v)  and  then  repeated. 

Figure  3  shows  (x,y)  scatter  plots  where  the  y' s  are  random  numbers  generated  by  the 
CRAY-1  computer  and  x,  —  (/  -  l/i)/  N. 


Figures  1,  2,  and  3  appear  to  show  that  the  Fibonacci  number  and  bit  reverse  methods 
produces  very  good  uniformity  and  that  the  random  number  generator  produces  considerable 
local  (and  probably  global)  bunching.  Let  us  define  a  goodness  number  D,  following  Wick 
(1979)  in  order  to  measure  the  uniformity  of  the  number  sets  just  produced. 


D  S  max 


N*y_ 

N 


-xy\ 


(3) 


for  all  x  and  y  which  satisfy 


0  <  x  ,  y  <  1  (4) 

and  Nxy  is  the  number  of  particles  in  the  sub-rectangle  (0,0),  (0,y),  (x,0),  (x,y),  i.e. , 

0  <  x,  <  x  and  0  y  (5) 

If  the  number  set  is  uniformly  distributed  and  N  approaches  infinity,  the  corresponding  number 
D  should  be  zero.  D  was  calculated  for  the  three  methods  mentioned  above  is  shown  in  Figure 
4  as  a  function  of  N.  The  Fibonacci  number  and  bit  reversed  methods  produce  essentially  the 
same  degrees  of  uniformity  with  D  =  3  /  N.  The  random  number  result  is  0  =  1/  -JJi . 
Hence,  Fibonacci  and  bit  reversed  are  smaller  by  roughly  3  / -JN . 

As  a  further  test,  we  used  the  cumulative  velocity  distribution  function  to  assign  velocities 
vx  to  a  set  of  beams  whose  envelope  is  /(vx),  and  then  used  the  three  different  number  sets  to 
scramble  the  particle  positions,  x.  For  a  given  Maxwellian  velocity  distribution  function,  the 
goodness  of  the  local  thermal  equilibria  were  examined,  by  calculating  the  average  of  the  first 
six  velocity  moments  over  a  width  Ax  Clothe  system  length.  We  then  compared  the  width  Ax 
at  which  the  average  local  velocity  moments  departed  by  10%  from  the  expected  values,  for 
each  of  the  three  different  number  sets.  We  found  that  for  a  large  particle  number  N(  ;>  103) 
both  the  Fibonacci  number  method  and  bit  reversed  method  produce  very  good  local  thermal 
equilibria.  Details  will  be  given  in  the  next  QPR. 
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42 

F.  ESI  Code 

No  special  progress  to  report. 


G.  EMI  Code 

No  special  progress  to  report. 


H.  EZOHAR  Code 

No  special  progress  to  report. 


I.  RINGHYBRID  Code 

No  special  progress  to  report. 
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SECTION  III 

SUMMARY  OF  REPORTS,  TALKS,  PUBLICATIONS,  VISITORS 

A.  Talks 

(a)  Three  talks  were  presented  at  the  APS  Division  of  Plasma  Physics  meeting  Nov.  10  -  14, 
1980,  San  Diego;  the  abstracts  were  given  in  the  previous  QPR. 

(b)  See  paper  by  Hamed  in  this  QPR,  presented  at  Compact  Torus  Meeting,  Dec.  2  -4,  1980, 
LASL. 

B.  Publications 

ERL  Memo.  No.  UCB/ERL  M80/40  by  Y-J  Chen  and  C.K.  Birdsall  was  accepted  by  Physics  of 
Fluids. 

ERL  Memo  No.  UCB/ERL  M80/20  by  Y-J  Chen  and  B.I.  Cohen  was  accepted  by  Physics  of 
Fluids. 


C.  Visitor 

M.J.  Gerver  spent  a  week  with  us  Nov.  1980. 


D.  Conference 

Workshop  on  Long  Time  Steps  in  Particle  Simulation,  Dec.  8,  9,  1980  was  organized  by  C.K. 
Birdsall  and  A.  Friedman  of  UCB  and  B.I.  Cohen  and  A.B.  Langdon  of  LLL.  The  meeting 
notice  and  schedule  follows: 
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WORKSHOP  CN  IONG  TIME  STEPS  IN  PARTICLE  SIMUIATION 


December  8  and  9  (Monday  and  Tuesday) ,  1980 
8:30  a.m.  to  5:30  p.m. 

Woodward  Roan 

Bechtel  Engineering  Center 

University  of  California,  Berkeley 

The  recent  work  of  Jacques  Denavit  and  Rod  Mason  on  implicit  integration 
schemes  and  of  Cohen  et  al. ,  on  orbit  averaging  has  inspired  this  meeting.  Many 
others  have  also  contributed  ideas  on  relaxing  the  numerical  stability  conditions 
on  the  time  step  in  simulations. 

The  meeting  format  allows  plenty  of  time  for  speakers  and  discussion.  The 
tentative  program  is  attached.  The  only  special  rule  for  the  workshop  is  that 
one  person  speak  at  a  time,  as  there  is  ample  time  for  all  viewpoints. 

The  objects  of  the  workshop,  in  no  special  order,  are: 

to  discuss  in  great  detail  Denavit  and  Mason's  implicit 
schemes  and  the  BAAL  implicit  algorithm 

TO  collect  ideas  on  u>  A»1  stability  and  accuracy 
requirements  and  kvAtpissues 

appropriateness  of  centered  vs.  damping  time  integration  algorithms 

to  address  possible  new  problems  associated  with  implicitness, 
e.g. ,  generalizations  to  multidimensions  and/or  EM  codes,  boundary 
conditions,  spatial  filtering 

to  choose  problem  areas  of  importance  for  applications 

to  understand  possible  enhanced  nosiness  in  implicit  methods, 
e.g.,  due  to  collection  of  higher  plasma  moments 

to  design  implicit  time  differencing  algorithms 

to  consider  melding  implicit  methods  with  time  filtering 
and  orbit-averaging  methods  in  order  to  reduce  noise 

lastly,  to  categorize  and  generalize  methods  that  do  not 
and  cannot  work 


The  meeting  rocm  is  a  few  steps  from  Cory  Hall  (northeast  comer  of  Campus) 
where  we  have  our  corputer  terminals  connected  to  JWFECC.  Henoe,  we  can  try 
things  out  at  any  time  during  the  workshop.  This  workshop  is  limited  to  the 
twenty  people  invited. 


Charles  K.  Birdsall 

Alex  Friedman 

EE3CS  Dept.,  Cory  Hall 

University  of  California,  Berkeley 

Berkeley,  CA  94720 


Bruce  Cohen 
A.  Bruce  Langdon 
Lawrence  Livermore 
National  Laboratory 


Schedule 
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Workshop  on  long  Tine  Steps  in  Particle  Simulation,  December  8,  9,  1980,  U. 


Monday,  8:30  a.m.  Welcoming  remarks  C.K.  Birdsall 

A.  Morning  Session  8:31  a.m.  A.  Bruce  Langdon,  Chairman 


Jacques  Denavit  "Time  Filtering  Particle  Simulations  with 

Northwestern  u  4»1" 

pe 

Rodney  Mason  "Implicit  Moment  Particle  Simulation  of  Plasmas" 

IASL 


Jack  Byers  "Long  Time  Step  Electrostatic  Quasineutral  Model 

T. IJi  for  Magnetized  Plasmas” 

Vinoent  Thomas  "Orbit  Averaging  in  Implicit  and  Explicit 
Bruce  Cohen  Electrostatic  Codes" 

U. C.  Berkeley 

Lunch  12:30  Faculty  Club,  Lewis  Roam 

B.  Afternoon  Session  1:30  p.m.  Jack  Byers,  Chairman 

Jerry  Brackbill  " implicit  Methods  Applied  to  2d,  EM  Codes" 

David  Forslund 
IASL 


A.  Bruce  Langdon  "BAAL  implicit  Particle  Algorithms" 
Alex  Friedman 
Bruce  Cohen 
ILL 

Alex  Friedman  "An  Experimental  BAAL  Code" 

Bruce  Oohen 
A.  Bruce  Langdon 

Bruce  Cohen  "Design  of  Time  Integration  Algorithms" 

A.  Bruce  Langdon 
Alex  Friedman 


C.  Berkeley 

Time  in 
Minutes 

60 

60 

15 

15 

60 

60 

60 

60 
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Tuesday,  8:30  a.m. 

C.  Morning  Session  8;  30  a.m.  A3am  Drcbot,  Chairman 


Time  in 
Minutes 


W.  W.  Lee  "Gyrokinetic  Particle  Pushing  Model"  30 

PPPL 

"Multiple  Scale  Simulation  Model" 

(Simulation  of  drift  wave  problems  in  a 
hcmogenious  plasna) 


Brerrian  Godfrey  "Iterative  Inversion  of  Implicit  Equations"  30 

MRC 

Viktor  Decyk  "Dawson  Long  Time  Step  Schemes"  30 

UCTA 


Bruce  Cohen  "Explicit  Magneto-Inductive  Hybrid  Code,  Orbit 

LLL  Averaging  and  Filtering" 


Robert  Freis  "2d  Orbit  Averaged  Simulation  of  a  Mirror  Machine"  30 

Bruce  Cohen 

rr.T. 


Lunch  12;  30  Faculty  Club,  Lewis  Room 

D.  Afternoon  Session  1:30  p.m.  C.  K.  Birdsall,  Organizer 

Rump  Sessions,  to  be  set  up  by  discussion  leaders  picked  during  Sessions  A,B,C, 
in  order  to  cover  topics  needing  more  discussion. 


End  by  5:30  p.m. 
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